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ABSTRACT 

We combine deep X-ray survey data from the Chandra observatory and the wide- 
area/shallow XMM-XXL field to estimate the AGN X-ray luminosity function in the 
redshift range z = 3 — 5. The sample consists of nearly 340 sources with either pho¬ 
tometric (212) or spectroscopic (128) redshift in the above range. The combination 
of deep and shallow survey fields also provides a luminosity baseline of three orders 
of magnitude, Lx {2 — lOkeV) « lO"*^ — 10‘*®ergs“^ at z > 3. We follow a Bayesian 
approach to determine the binned AGN space density and explore their evolution in 
a model-independent way. Our methodology properly accounts for Poisson errors in 
the determination of X-ray fluxes and uncertainties in photometric redshift estimates. 
We demonstrate that the latter is essential for unbiased measurement of space den¬ 
sities. We find that the AGN X-ray luminosity function evolves strongly between the 
redshift intervals z = 3 — 4 and z = 4 — 5. There is also suggestive evidence that 
the amplitude of this evolution is luminosity dependent. The space density of AGN 
with Lx{2 — lOkeV) < lO^^ergs”^ drops by a factor of 5 between the redshift inter¬ 
vals above, while the evolution of brighter AGN appears to be milder. Comparison of 
our X-ray luminosity function with that of UV/optical selected QSOs at similar red- 
shifts shows broad agreement at bright luminosities. Lx (2 — 10 keV) > lO"*® erg s“^. At 
fainter luminosities X-ray surveys measure higher AGN space densities. The faint-end 
slope of UV/optical luminosity functions however, is steeper than for X-ray selected 
AGN. This implies that the type-I AGN fraction increases with decreasing luminosity 
at z > 3, opposite to trends established at lower redshift. We also assess the signifi¬ 
cance of AGN in keeping the hydrogen ionised at high redshift. Our X-ray luminosity 
function yields ionising photon rate densities that are insufficient to keep the Universe 
ionised at redshift z > 4. A source of uncertainty in this calculation is the escape 
fraction of UV photons for X-ray selected AGN. 

Key words: galaxies: active - galaxies: Seyferts - X-rays: diffuse background - 
quasars: general 
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1 INTRODUCTION 

In recent years observations have established that super- 
massive black holes (SMBHs) are nearly ubiquitous in local 
spheroids (Magorrian et al. [1998 Kormendy & Ho 20131. 
These relic black holes are believed to have grown their 
masses at earlier times mostly via accretion of material from 
larger scales (e.g. Soltan|[l982 Marconi et al.p004 l. Ques¬ 
tions that remain open are when during the lifetime of the 
Universe these events occurred and under what physical con¬ 
ditions black holes grow their masses. Moreover, observa¬ 
tions show that in the local Universe correlations exist be¬ 
tween the mass of SMBHs and the properties of the stellar 
component of the bulges in which they reside, such as ve¬ 
locity dispersion (e.g [Ferrarese fc Merritt||2000| |Gebhardt| 


at faint accretion luminosities. Observations at X-ray wave¬ 


lengths can mitigate these issues (e.g Brandt & Alexan 
der|2015 K X-ray photons, particularly at rest-frame energies 


et aL|2000 

Giiltekin et al.|2009 Graham et aL|2011 1, lumi- 

nosity (e.g 

McLure & Dunlop||2002| Marconi & Hunt||2003 

Giiltekin et al.|2009|), dynamical mass (e.g. Magorrian et al. 

1998 

[Marconi & Hunt||2003| [Haring & Rix 2004 Graham 


Savorgnan et al.|2013[). Such correlations suggest a link be¬ 


tween the growth of black holes and the formation of their 
host galaxies, although the exact nature of such an interplay 
is still not well understood. Processes that can establish such 
correlations include a common gas reservoir that both feeds 
the central black hole and forms stars on larger scales, out¬ 
flows related to the energy output from the AGN itself that 
affect the Inter-Stellar Medium and regulate the formation 


> 2 keV, can penetrate nearly unaffected large columns of 
intervening gas and dust clouds (Nh ^ 10^^ cm“^), thereby 
providing samples least affected by obscuration biases. 
Moreover, the X-ray emission associated with stellar pro¬ 
cesses is typically 2 orders of magnitude fainter than the 
AGN radiative output and therefore contamination or dilu¬ 
tion effects by the host galaxy are negligible at X-rays over 
a wide baseline of accretion luminosities. X-ray surveys also 
benefit from a well defined selection function that is rela¬ 
tively easy to quantify and account for in the analysis. The 
disadvantage of X-ray selection is that the detected AGN 
are often optically faint and therefore spectroscopic follow¬ 
up observations are often expensive. Nevertheless, intensive 
multiwavelength campaigns in recent years substantially in¬ 
creased the number of X-ray survey fields with sufficient 
quality ancillary data for reliable X-ray source identification 
and redshift measurements using either spectroscopy or pho¬ 
tometric methods. Early results by |Cowie et al. (2003 \ on 
the redshift evolution of the X-ray AGN space density and 


Ueda et al. (20031 on the hard band (2-10 keV) X-ray lumi¬ 


nosity function and obscuration distribution of AGN have 
been expanded recently both in terms of data and analysis 


of stars (Silk & Rees 

1998 Fabian 199^ 

King 2003 2005 

|Di Matteo et al.||2005 

[Groton et al.||2006 

), or the merging 

history of galaxies and their SMBHs (e.g. Jahnke & Maccio 


2011 ). 


methodology ( 

Yencho et al. 

2009 Ebrero et al. 

2009 

Aird 

et al. 

|2010[ [Burlon et al.||2011[ [Ueda et al.||2014[ [Buchner 

et al. 

|2015| [Aird et al.||2015 

Miyaji et al.||2015 

. Although 

important details on the shape of the X-ray luminosity func¬ 
tion and the obscuration distribution of AGN are still de- 

bated (Ueda et al. 2014 Buchner et al. 2015 

Aird et al. 


One approach for improving our understanding of the 
formation of SMBHs as a function of cosmic time and their 
relation to their host galaxies is population studies of Ac¬ 
tive Galactic Nuclei (AGN), which signpost accretion events 
onto SMBHs. This requires a census of the AGN population 
across redshift to constrain for example, the accretion his¬ 
tory of the Universe or study the incidence of active black 


luminosity function is reasonably well constrained at least 
to j* « 3. The initial increase of the AGN X-ray luminosity 
density from z = 0 to z ~ 1.5 is followed by a broad plateau 
up to z « 2.5 — 3 and a decline at higher redshift. However, 
the amplitude of the X-ray AGN evolution at 2^3 is still 
not well constrained. Early studies suggested a moderate de¬ 


holes among galaxies. In that respect, the AGN luminosity 

Ebrero et al.|2009 Aird et al. 

2010), contrary to claims for a 

function, i.e. their comoving space density as a function of 

rapid drop (Brusa et al.|2009 

Civano et al.|2011 Vito et al. 

redshift and accretion luminosity, is one of the fundamental 

|2013||Kalfountzou et al.|2014|) that can be parametrised by 

quantities that characterise the demographics of active black 

an exponential law in redshift (Gilli et al. 2007) similar to the 

holes. The cosmic evolution of AGN leaves imprints on the 

optical QSO space density evolution (Schmidt et al. 1995 


shape and overall normalisation of the luminosity function. 
The total mass density locked into black holes at different 
epochs can be inferred by direct integration of the AGN lu¬ 
minosity function, under assumptions about the radiative 
efficiency of the accretion process and after applying appro- 


small X-ray AGN sample sizes at 2 > 3. For example, there 
are 209 and 141 X-ray AGN with spectroscopic or photomet¬ 
ric redshifts above 2 = 3 in the most recent compilations of 


priate bolometric luminosity corrections (e.g Marconi et al. 
2004 Aird et al.|2010 Ueda et al.|2014 1. The space density 


Kalfountzou et al. (2014) and Vito et al. (2014) respectively. 


of AGN split by host galaxy properties, such as stellar mass, 
morphology or level of star-formation, provides clues on the 
interplay between black hole accretion and galaxy evolution 
( Georgakakis et al.||2009 2011 Aird et n(]|2012 Bongiorno 


These numbers should be compared with sample sizes of few 
thousands AGN at 2 < 3 for the most recent X-ray lumi¬ 
nosity function studies (e.g Ueda et al.|2014 Buchner et al. 


et al.|20T2 Georgakakis et al.||2014 ). 


Selection at UV/optical wavelengths (e.g. Richards 
et al.|[2009 Ross et al.||2012 ) currently provides the largest 
spectroscopic AGN samples for luminosity function calcu- 


2015 Aird et al.|2015 Miyaji et al.|2015[ 

Better constraints on the form and amplitude of the 
evolution of X-ray AGN at 2 > 3 will have implications 
for the contribution of this population to the UV photon 
field density that is needed to keep the Universe ionised 


lations (e.g. Ross et al. 2013). The downside is that the 
UV/optical continuum of AGN is sensitive to dust extinc¬ 
tion along the line-of-sight and dilution by the host galaxy 


at high redshift. Haardt & Madau (2012) used the Ueda 


et al. (2003) X-ray luminosity function to predict a moder¬ 


ate contribution of AGN to the hydrogen ionising radiation 
field at 2 > 3, in broad agreement with constraints derived 
from UV/optical selected QSO luminosity functions (e.g. 
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The AGN XLF at z = 3 — 5 3 


Fontanot et al.||2007| [Masters et af.||2012| [McGreer et af.| 

2013 

I and other X-ray AGN studies ( 

Barger et al. 

2003a 

Grissom et al. 20141. There are also 

claims however, that 
n to the photoionisa- 
2012 Glikman et al. 

AGN provide an important contributio 
tion rate at high redshift (Fiore et al.| 


2011 Giallongo et al. 20151. This discrepancy emphasises 


the need for further work to improve measurements of the 
AGN space density at high redshift and to better understand 
their role in the re-ionisation of the Universe. 

In this paper we combine deep Chandra and wide- 
area/shallow XMM-Newton survey fields to compile one of 
the largest samples of X-ray selected AGN at 2 = 3 — 5 
to date. A Bayesian methodology is developed to correctly 
account for photometric redshift uncertainties and to deter¬ 
mine in a non-parametric way the AGN comoving space den¬ 
sity in the redshift intervals 2 = 3—4 and 2 = 4—5 and over 3 
decades in X-ray luminosity [log Lx (2 —10 keV) « 43 — 46 in 
erg s“ . Although parametric models are also fit to the data, 
we emphasise the importance of non-parametric estimates to 
determine in a model-independent way the shape and over¬ 
all evolution of the X-ray luminosity function. Throughout 
this paper we adopt Ho = 70kms“^ Mpc“^, Hm = 0.3 and 
flA = 0.7. 


2 DATA 

For the determination of the X-ray luminosity function in 
the redshift interval 2 = 3 — 5 we combine Chandra and 
XMM-Newton X-ray surveys with different characteristics in 
terms of area coverage and X-ray depth. These are the 4 Ms 


et al. 20131, the 2 Ms Chandra Deep Field North (CDFN; 

Alexander et al. 2003 Rangel et al. 2013|), the Extended 

Groth Strip International Survey field 

(AEGIS, [Davis et al.| 

2007 Laird et al.||2009 Nandra et al. 

20151, the Extended 

Chandra Deep Field South (ECDFS; 

Lehmer et al. 20051, 


and the equatorial field of the XMM-XXL survey. 


2.1 Chandra survey fields 

The Chandra observations of the CDFS, CDFN, AEGIS, 
ECDFS and C-COSMOS were analysed in a homogeneous 
way by applying the reduction and source detection method¬ 


ology described by Laird et al. (20091. Specific details on 


the analysis of the 4 Ms CDFS and the 2 Ms CDFN fields 


are presented by Rangel et al. (20131. The Chandra sur¬ 
vey of the AEGIS field has two tiles. The wide and shallow 
one (AEGIS-W) consists of 8 Chandra pointings of 200 ks 
each. These data are described by Laird et al. (20091. The 


deep survey of the AEGIS field (AEGIS-XD) increased to 
a total of 800 ks the exposure time of the central regions of 
the AEGIS-W. The additional data overlap with the central 
3 of the original 8 Chandra pointings observed as part of 
the AEGIS-W. The AEGIS-XD survey data reduction and 
source catalogue generation are described by|Nandra et al.| 
( 201 ^. 


The X-ray sources used in this paper are detected in the 
0.5-7 keV (full) spectral band with Poisson false detection 
threshold < 4 x 10“®. The count rates in the 0.5-7 keV band 
are converted to fluxes in the 0.5-10 keV band assuming a 


power-law spectral index with F = 1.9. This is steeper than 
the F = 1.4 adopted for the X-ray flux estimation in the 
published catalogues of the CDFS, CDFN and AEGIS fields 
( Laird et al.||2009 Rangel et al.|[20l3 Nandra et al.|[20T5 |. 
The choice of the F = 1.9 is motivated by the fact that 
at high redshift, 2 : >3, the observer-frame 0.5-7keV band 
corresponds to harder rest-frame energies, which are least 
affected by obscuration. A diagnostic of the spectral shape 
of X-ray sources is their hardness ratio defined as HR=(H- 
S)/(H-|-S), where S, H are the observed count rates in the 
0.5-2 and 2-7 keV spectral bands, respectively. For sources 
in the range 2 = 3 — 5 we find a median hardness ratio HR ~ 
—0.2. This is consistent with a power-law X-ray spectrum 
with F « 1.8, i.e. similar to the mean spectral index of the 
intrinsic AGN spectra (F « 1.9, Nandra fc Pounds||1994 1. 
We therefore choose to fix F = 1.9 for the determination of 
fluxes. We note however, that the choice of F (1.4 vs 1.9) 
has a small impact on the results presented in this paper. 

Sensitivity curves, which measure the total survey area 
that is sensitive to sources of a particular flux are calculated 


following methods described in Georgakakis et al. (20081. 


The overlap between the ECDFS and the 4 Ms CDFS or the 
AEGIS-W and the AEGIS-XD is accounted for by defining 
independent spatial regions for each survey. Spatial masks 
that describe both the boundaries of the optical/infrared 
imaging of each field and regions of poor photometry because 
of bright stars ( Aird et al.|2015 1 are also taken into account 
in the X-ray sensitivity calculations. The sensitivity curves 
in the 0.5-10 keV band are presented in Figure Table 
presents the number of X-ray sources in each field. The same 
spatial masks used for the construction of sensitivity maps 
are used to filter the X-ray source catalogue. 

The optical identification of the X-ray sources in the 
CDFN, AEGIS-XD. AEGIS-W, ECDFS and C-COSMOS 


fields are based on the Likelihood Ratio method (Suther- 
land fc Saunders|1992 1 as implemented in Aird et al. (20151. 


The multiwavelength associations of the 4 Ms CDFS X-ray 


sources are presented by Hsu et al. (20141. They apply 


a Bayesian methodology, based on the work of |Budavari] 
& Szalay (20081, to different catalogues available in that 


field including the CANDELS/iL-band selected photome¬ 
try presented by |Guo et al.| ( |2013p, the T aiwan ECDFS 
Near-Infrared Survey (TENIS; |Hsieh et al. 20121 and the 
MUSYC/BUR-selected catalog of Cardamone et al.| ( |2010 |. 
The number of X-ray sources with secure optical or infrared 
counterparts in each field are presented in Table 

Extensive spectroscopic campaigns have been carried 
out in the fields of choice. For the CDFN, ECDFS and 
AEGIS-W we use the compilation of spectroscopic redshifts 


presented by Aird et al. (20151. In the 4 Ms CDFS we use 


the spectroscopic redshifts compiled by Hsu et al. (20141. 


In the case of AEGIS-XD we use the spectroscopic red¬ 


shift catalogue presented by (Nandra et al. 20151. Red¬ 


shifts in the C-COSMOS are from the public releases of 


the VIMOS/zCOSMOS bright project (Lilly et al. 20091 


and the Magellan/IMACS observation campaigns (Trump 


et al. 20091, as well as the compilation of redshifts for X- 
ray sources presented by |Civano et al. (20121. The spectro¬ 
scopic redshifts used in this paper have quality flags in the 
published catalogues from which they were retrieved that 
indicate a probability better than « 95% of being correct. 

For X-ray sources without spectroscopy, photometric 
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4 Georgakakis et al. 


Table 1. Number of X-ray sources in the full-band selected sample 


field 

solid angle 

Number of 

Number of 

photometric redshift 


spectroscopic redshift 



(arcmin^) 

X-ray sources 

optical/infrared IDs 

full sample 

3 < 2 < 4 4 < 

; 2 < 5 

full sample 

3<2<4 4<2<5 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 

4Ms CDFS 

271.4 

422 

418 

151 

16 

3 

268 

9 

1 

CDFN 

412.9 

453 

436 

171 

11 

5 

265 

9 

2 

ECDFS 

643.0 

407 

399 

267 

30 

11 

134 

3 

1 

AEGIS-XD 

934.6 

818 

804 

478 

37 

13 

326 

16 

0 

AEGIS-W 

724.3 

415 

413 

267 

26 

8 

146 

2 

1 

COSMOS 

3258.0 

1435 

1393 

635 

38 

14 

758 

21 

4 

XMM-XXL 

64728.0 

7493 

3798 

- 

- 

- 

2553 

54 

5 

Total 

70972.2 

- 

- 

- 

158 

54 

- 

114 

14 


(1) Name of the X-ray survey fields used in this work. The last row lists the total number of sources for the combined fields. For that 
last row only the columns that correspond to the total number of sources with photometric or spectroscopic redshift in the intervals 
3 < 2 < 4 and 4 < 2 < 5 are listed. This is because these are the relevant numbers to the analysis presented in this paper. (2) Solid 
angle of each sample in square arcminutes after excluding regions of poor/no photometry. (3) Total number of X-ray sources detected 
in the 0.5-7 keV energy band in the case of the Chandra, or the 0.5-8 keV band for the XMM-XXL sample. (4) Number of full-band 
detected X-ray sources with optical/infrared associations. (5) Number of sources with photometric redshift estimates in each survey field. 
Photometric redshifts for the XMM-XXL are not available. (6) Number of sources with photometric redshifts in the interval 2 = 3 — 4. 
This the sum of the photometric redshift Probability Distribution Functions rounded to the nearest integer. No photometric redshift 
estimates are available for the XMM-XXL sample. (7) Same as column (6) but for the redshift interval 2 = 4 — 5. (8) Number of sources 
with spectroscopic redshift estimates in each survey field. (9) Number of sources with spectroscopic redshifts in the interval 2 = 3 — 4. 
(10) Number of sources with spectroscopic redshifts in the interval 2 = 4 — 5. 



(0.5-10keV) (erg s“* cm“^) 



(0.5-10keV) (erg s“'cm“®) 


Figure 1. X-ray Sensitivity curves for the combined Chandra surveys used in the analysis (left panel) and the XMM-XXL survey (right 
panel). 


redshifts are estimated using the multiwavelength photo¬ 
metric catalogues available for each survey held. The pho¬ 
tometric redshifts of the X-ray sources in the 4 Ms CDFS, 
the AEGIS-XD and the COSMOS helds are determined fol¬ 


lowing the methodology described by Salvato et al. (2009 


2011). Specihc details can be found in Hsu et al. (2014; 4 Ms 


CDFS), Nandra et al. (2015; AEGIS-XD) and Salvato et al. 
(2011; COSMOS held). The estimated rms scatter of the X- 
ray AGN photometric redshifts is crAz/(i+z) = 0.016, 0.014 


and 0.04 for the C-COSMOS, 4 Ms CDFS and AEGIS-XD 
samples, respectively. The corresponding outlier fraction, 
dehned as /S.z/{l+z) > 0.15, is about 5-6% in all three helds. 
In the case of ECDFS, AEGIS-W and CDFN we use the 


photometric redshifts estimated by |Aird et al.| | |2015[ ) . For 
these helds o'az/(i+z) = 0.06 and the outlier fraction is about 
15%. The latter value is larger than in the C-COSMOS, 4 Ms 
CDFS and AEGIS-XD helds. This is related to diherences 
in the methodology of estimating photometric redshifts and 
ultimately to the choice of template SEDs used in the cal¬ 
culation. Nevertheless, the photometric redshifts estimated 
by Aird et al. (20151 are assigned appropriately larger un¬ 
certainties, approximated by the corresponding Probability 
Distribution Functions (PDZ), that rehect the higher out¬ 
lier fraction. Figure plots spectroscopic vs photometric 
redshifts for the sample used in this paper and illustrates 
the overall quality of the photometric redshifts estimates. 
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Figure 2. Spectrocopic vs photometric redshift measurements for 
the sample of X-ray selected sources used in this paper. The data 
points correspond to the median value of the PDZ. The errorbars 
correspond to the 90% confidence interval around the median. 


Figure 3. Examples of photometric redshift Probability Distri¬ 
bution Functions for X-ray selected AGN in the AEGIS field. 
The different colours correspond to different sources. There is a 
variety of PDZ shapes in the sample, including uni-modal and 
relatively narrow (blue solid), unimodal but broad (red dashed) 
and multi-modal (black dot-dashed). 


In the X-ray luminosity function calculations we use the 
full photometric redshift PDZ.These are typically unimodal 
but at increasing redshift they broaden and secondary peaks 


may also appear. Also the Aird et al. (20151 PDZs are typ¬ 
ically broader than those in the C-COSMOS, 4 Ms CDFS 
and AEGIS-XD fields. Examples of PDZs used in this paper 
are shown in Figure]^ Tablepresents the number of pho¬ 
tometric and spectroscopic redshifts in each field, both total 
and in the intervals 2 = 3 — 4 and 4 — 5. Sources without 
optical identifications and hence, without redshift estimates, 
are a minority in the Chandra surveys sample, 107 in total. 
These sources can be either moderate redshift (2 « 1 — 3) 
AGN with red SEDs because of e.g. obscuration and/or old 
stellar populations ([Koekemoer et al.||2004 Schaerer et al. 


|2007[ [Rodighiero et al.||2007[ |Del Moro et ai.|2009p , or high 

redshift systems (2 >3). Since the redshift distribution of 
these sources is not known they are assigned a flat PDZ in 
the redshift interval 2 = 1 — 6 and zero at other redshifts. 
Fixing the lower limit of the redshift range above to a value 
between 2 = 0 and 2 = 2 does not change the results and 
conclusions. 

Additionally we have tested that the differences in the 
accuracy and outlier fraction of the photometric redshifts 
in the different fields used in this work do not affect the 
final results. For the COSMOS, AEGIS-XD and 4 Ms CDFS 
we can substitute the photometric redshift PDZs adopted 
in this paper (based on the methods of Salvato et al. 2009, 
2011) with those estimated following the methodoloy of Aird 
et al. (2015). For these fields we can therefore estimate the 
X-ray luminosity function at 2 = 3 — 5 (see next sections) 
using two different sets of photometric redshifts, those of 
Aird et al. (2015) and those determined following Salvato et 
al. (2009, 2011). We find no systematic differences between 
the X-ray luminosity functions at 2 = 3 — 5 estimated using 
the two independent photometric redshift catalogues. 


2.2 XMM-XXL survey data 

The Chandra survey fields are complemented by the wide- 
area and shallow XMM-XXL survey (PI: Pierre). The XMM- 
Newton observations of this programme cover a total of 
about 50 deg^ of sky split into two nearly equal area fields. 
In this paper we use the equatorial region of the XMM-XXL, 
which overlaps with the Canada-France-Hawaii Legacy Sur¬ 
vey (CFHTLS) W1 field and covers an area of about 25 deg^. 
The approximate centre of this field lies at right ascension 
a = 10 : 22 : 42.20 and declination S = —04 : 52 : 51.03. 

The reduction of the XMM data, the construction of 
the X-ray catalogue and the association of the X-ray sources 


with optical counterparts are described by Liu et al.| (2015 \ 


based on the methods presented in [Georgakakis fc Nandr^ 
| |2011[ ). In brief, the X-ray data reduction is carried out us¬ 
ing the XMM Science Analysis System (SAS) version 12. We 
analyse XMM-Newton observations related to the XMM- 
XXL programme that were made public prior to 23 January 
2012. XMM-XXL data observed after that date are not in¬ 
cluded in the analysis. As a result our final catalogue of 
the equatorial XMM-XXL field misses about 5deg^ worth 
of X-ray coverage. The epchain and emchain tasks of SAS 
are employed to produce event files for the EPIC (Euro¬ 
pean Photon Imaging Camera; [Striider et al.|[2001| | Turner | 
et al.|200T I PN and MOS detectors respectively. Flaring pe¬ 
riods resulting in elevated EPIC background are identified 
and excluded using a methodology similar to that described 


by Nandra et al. (2007). We use X-ray sources detected 


in the 0.5 — 8 keV spectral band with Poisson false detec¬ 
tion probability of < 4 x 10“®. The final sample consists 
of 7493 unique sources detected in 0.5-8 keV spectral band. 
The fluxes listed in the final source catalogue are in the 0.5- 
10 keV band assuming a power-law spectral energy distribu¬ 
tion with P = 1.4. These X-ray sources are matched to the 
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SDSS-DR8 photometric catalogue ( [Aihara et aI.|20lT I using 
the Maximum-Likelihood method (Sutherland & Saunders 
1992). We assign counterparts to 3798 sources with Likeli¬ 
hood Ratio LR > 1.5. At that cut the spurious identification 
rate is about 6% and the total number of 0.5-8 keV detected 
sources with optical counterparts is 3798 (see Table [^. 

Redshifts for the XMM-XXL X-ray sources are from 
several follow-up spectroscopic campaigns. The XMM- 
XXL field overlaps with the SDSS-III Baryon Oscillation 


Spectroscopic Survey (BOSS; Dawson et al.| 20131 pro¬ 
gramme, which provides spectroscopy for UV/optically se¬ 
lected broad-line QSOs and luminous red galaxies. |Stalin| 
et al. (20101 presented spectroscopy for X-ray sources se¬ 
lected in the original XMM-LSS survey ( Clerc et al.||2014 |, 
which is part of the equatorial XMM-XXL survey field. 
Most of the redshifts however, are from a total of five spe¬ 
cial SDSS plates dedicated to follow-up spectroscopy of X- 
ray sources as part of the Ancillary Programs of SDSS-III. 
The overlap between those plates and the XMM-XXL sur¬ 
vey region is 17.98 deg“^. Targets were selected to have 
fx(0.5— 10keV,r = 1.4) > 10“^'*ergs“^ cm“^ and 15 < 
r < 22.5, where r is either the SDSS PSF magnitude in the 
case of optically unresolved sources (SDSS type=6) or the 
SDSS model magnitude for resolved sources. Specific details 
on these spectroscopic observations, including spectral clas¬ 
sification, visual inspection and redshift quality assessment 


are presented by Menzel et al. (20151. The total number of 


XMM-XXL sources with secure spectroscopic measurements 
are presented in Table Also shown in this table are the 
number of sources with spectroscopic redshifts in the inter¬ 
val 2 = 3 — 5. 


Although when constructing the X-ray source catalogue 
of the XMM-XXL field fluxes are estimated for a power-law 
X-ray spectrum with spectral index P = 1.4, in the rest of 
the analysis we adopt P = 1.9 for the calculation of fluxes, 
luminosities and sensitivity maps. This is because at the 
depth of the XMM-XXL survey AGN at 2 > 3 are powerful 
QSOs with Lx(2 — lOkeV) > 10^"* ergs“^. The fraction of 
obscured AGN among such luminous sources is a decreasing 
function of luminosi ty (|Ueda et al.|2003| |Akylas et aE 2006| 


Merloni et al.|[^014 Ueda et al.||2014' Buchner et al. 2015| 

and therefore a spectral index of P = 1.9, which represents 
the intrinsic unobscured power-law X-ray spectrum of both 



Figure 4. Monochromatic 2keV X-ray luminosity, L^(2keY), 
plotted as a function of monochromatic 2500 A UV luminosity, 
Li/( 2500A). The data points (cross or dots) are XMM-XXL X- 
ray selected broad line QSOs with secure spectroscopic redshifts 
in the interval 2 = 1 — 5. Sources at 2 > 3 are highlighted with dif¬ 
ferent symbols (crosses). The red dashed line is the bisector best- 
fit Li, (2 keV) — Li, (2500 A) relation determined by Lusso et al. 
[2010; logLakeV = 0.760 logLi/(2500 A) + 3.508]. For the XMM- 
XXL QSOs the X-ray luminosity density at 2 keV is estimated 
from the 0.5-2 keV flux assuming a power-law X-ray spectrum 
with index F = 1.9. The 2500 A monochromatic luminosity is 
determined from the SDSS photometry. For the k-corrections we 
adopt the simulated QSO templates of |McGreer et ar] ( |2013| |. For 
a QSO with redshift 2 the SDSS photometric bands with effective 
wavelengths that bracket the wavelength 2500 X (I-I- 2 ) A are iden¬ 
tified. The mean model QSO SED at that redshift is then scaled 
to the observed SDSS optical magnitudes in those two bands. The 
monochromatic luminosity at 2500 A is then estimated from the 
scaled model SED. At redshifts 2 > 2.7 the rest-frame 2500 A lies 
beyond the effective wavelength of the SDSS 2 -band (9134A). 
For these sources the flux density at 2500 X (1 + 2 ) A is an ex¬ 
trapolation using the model SED. The results do not change if 
we simply linearly interpolate between the observed flux densi¬ 
ties of the SDSS bands that bracket the rest-frame 2500 A. This 
model-independent approach however, does not allow extrapola¬ 
tion beyond 2 > 2.7. 


local Seyferts (e.g Nandra & Pounds] 1994 

and luminous 

high-redshift QSOs (e.g. Vignali et al.|2005 

Shemmer et al. 


The 0.5-10 keV fluxes estimated for P = 1.9 are about 35% 
fainter than those for P = 1.4. The XMM-XXL sensitiv¬ 
ity cnrve in the 0.5-10 keV band is shown i n Figure [T] and 
is estimated following methods described in [Georgakakis fc| 


Nandra (20111. In the calculation of the sensitivity curve 


we only consider the overlapping area between the SDSS- 
III Ancillary Programs spectroscopic plates used to target 
X-ray sources and the XMM-XXL survey region. We also 
take into account the flux limit for follow-up spectroscopy 


3 METHODOLOGY 

3.1 X-ray Luminosity function estimation: 
Chandra survey fields 

A Bayesian approach is adopted for the determination of the 
X-ray luminosity function. The X-ray sources detected in a 
survey are essentially Poisson realisations of a parent sample 
and therefore the likelihood can be written as the product 


/x(0.5 — 10keV,r = 1.4) > 10“^"^ erg s“^ cm“^. This limit 

of the Poisson probabilities of individual sources. Following 

appears as a smooth drop in area rather than a sharp cut in 

the works of Marshall et al. ( 

1983), Loredo (2004), Aird 

Figure [^because the Poisson probability of measuring a flux 

et al. (20101 and Buchner et al. 

(20151 the likelihood can be 

above this limit is used to determine the sensitivity curve. 

written as 
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logL^(2500A) - <logL^ (2500A)> 


Figure 5. Histogram of the difference A log L^, (2500 A) be¬ 
tween the monochromatic 2500 A luminosity of broad-line QSOs 
in the XMM-XXL (blue points of Figure and the mean < 
logLiy(2500 A) > value predicted by the bisector best-fit relation 
of Lusso et al. (2010; red dashed line of Figure for a given 
logLt/(2keV). The red dashed line shows the best-fit Gaussian 
distribution. The mean value of that distribution is consistent 
with zero and the standard deviation is cr = 0.4. 



Figure 6. Fraction of XMM-XXL X-ray sources with success¬ 
ful redshift measurement as a function of the r-band magnitude 
of the optical counterpart. The fraction is defined as the ratio 
between the number of potential targets [/x(0.5 — lOkeV, F = 
1.4) > 10“^^ erg s“^ cm~^ and secure optical counterparts with 
r < 22.5 mag] and X-ray sources with successful redshift measure¬ 
ments and fluxes/magnitudes within the above cuts. 



Figure 7. Impact of different levels of obscuration on the 0.5- 
lOkeV flux of a source as a function of redshift. The vertical 
axis is the ratio of the flux of an AGN obscured by column den¬ 
sity Nfj relative to the flux of the same source in the case of 
zero obscuratrion. The red curves correspond to different levels 
of obscuration in the range logTV^ = 22.5 — 24.5 (cm“^ units). 
The vertical dashed lines mark the redshift interval of interest, 
2 = 3 — 5. For fixed obscuration and intrinsic luminosity the 
flux of higher redshift AGN is less affected. The differential flux 
suppression between 2 = 3 and 2 = 5 is small, < 10%. For the 
calculation of X-ray k-corrections we adopt the model X-ray spec¬ 
tral of |Brightman & Nandr^ l |2011a1 l. These are based on Monte 
Garlo simulations of an illuminating source at the centre of a 
sphere with constant density and a conical region (apex at the 
centre of the sphere) cut-off to approximate a toroidal geometry. 
These simulations take into account both Gompton scattering and 
photoelectric absorption of the X-ray photons by the obscuring 
medium. We adopt F = 1.9 for the intrinsic AGN spectrum, an 
opening angle of the conical region of 60 deg and a viewing angle 
of 45 deg, i.e. a line-of-sight intersecting the obscuring material. 


C{d, I e) = e~^x 

^ r c\v ( 1 ) 

P[ / d log Lx -^dz p{di\ Lx,z) (l){L:^,z \ 6), 

where dV/dz is the comoving volume per solid angle at red¬ 
shift z, di signifies the dataset and 6 represents the param¬ 
eters of the luminosity function model, <f){Ly^,z \ 9), that 
are to be estimated. The multiplication is over all sources, 
N, and the integration is over redshift and X-ray luminos¬ 
ity. The quantity p{di\ Lx,z) is the probability of a par¬ 
ticular source having redshift z and X-ray luminosity Lx. 
This captures uncertainties in the determination of both red- 
shifts (e.g. photometric redshifts measurements) and X-ray 
fluxes because of Poisson statistics and the Eddington bias. 
In equation A is the expected number of detected sources 
in a survey for a particular set of model parameters 6 

\ = j dlogLx^d2^(Lx,«) <^(Lx,2 I 0). (2) 
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where, A[Ly^,z) is the sensitivity curve that quantifies the 
survey area over which a source with X-ray luminosity Lx 
and redshift z (and hence flux fx) can be detected. Note 
that the selection function term, A{Lx,z), is not included 
within the integral of equation 1. The reader is referred to 
an extensive discussion in Loredo (2004) on that point. 

The goal of this paper is the estimation of the X-ray 
luminosity function in the redshift interval 2 = 3 — 5. How¬ 
ever, defining a sample of X-ray AGN in a relatively nar¬ 
row redshift range is not straightforward. This is because 
the redshifts of many sources are determined by photomet¬ 
ric methods and therefore have uncertainties, which are not 
negligible compared to the size of the redshift interval. It 
may happen for example, that the errors of the photometric 
redshift of a particular source straddle one (or both) bound¬ 
aries of the redshift range of interest (see Fig. [^. 

We deal with this difficulty by simply using all sources 
in the X-ray sample and splitting the luminosity function 
model into two terms 


(piLx, z\9) = z G [zi,Z 2 ] | 9i)+ 

<I>2{Lx,z ^ [ 21 , 22 ] I 92 ). 


(3) 


The first term (j>i refers to the luminosity function within the 
redshift interval of interest, 2 = 21 — 22 , and has its own set 
of parameters, 9i. The second term, (j) 2 , corresponds to the 
X-ray AGN space density outside that redshift range and 
has a different set of parameters, 02 , which are treated as 
nuisance parameters. The advantage of this approach is that 
it allows an estimate of the AGN space density at 2 = 3 — 5 
nearly independent of the shape and form of the evolution of 
the X-ray luminosity function at lower redshift. In any X-ray 
selected sample the bulk of the AGN population lies at low 
redshift 2 < 3. This may introduce systematics in the deter¬ 
mination of the AGN space density at 2 > 3, if the same evo¬ 
lutionary law is fit to the data across all redshifts. In this case 
it is possible that the determination of the model parame¬ 
ters is dominated by the regions of the redshift/luminosity 
parameter space with the most data. Our approach min¬ 
imises the impact of this potential source of bias. In our 
analysis the (f >2 term is modelled as a step function, i.e. the 
sum of constants that correspond to the AGN space densi¬ 
ties at different luminosity and redshift bins. The values of 
these constants are determined by the data via equation 
For the calculation of equation assumptions need to be 
made on the redshift errors. For sources with photometric 
redshift determinations we adopt the corresponding redshift 
Probability Distribution Function (PDZ) as a measure of 
the uncertainty. Spectroscopic redshifts in the sample have 
reliabilities > 95% and therefore their corresponding PDZs 
are assumed to be delta functions at the spectroscopic red¬ 
shift of the source. Sources without optical counterparts are 
assigned a flat PDZ in the redshift range 2 = 1 — 6 (see 


Section 2.11. 


Poisson statistics are used to determine the flux distri¬ 
bution that is consistent with the extracted source and back¬ 
ground counts. A power-law X-ray spectrum with F = 1.9 is 
adopted in this calculation. The flux distribution in the 0.5- 
10 keV band is then convolved with the PDZ to estimate the 
luminosity distribution of each source at rest-frame energies 
2-10 keV. The relevant k-corrections also assume a power- 
law X-ray spectrum with P = 1.9. The two-dimensional 


probability distribution in Lx and 2 is the term p{di\ Lx, z) 


|et al.|l992[ | to evaluate the integral of equationFor each 
source we draw Lx and 2 samples based on the PDZ and 
Poisson X-ray counts distribution of that source. The lu¬ 
minosity function is then evaluated for each sample point, 
Lx,z. The integral of equation]^ is simply the average lu¬ 
minosity function of the sample. 


of equationlH In practice we use importance sampling (Press 


3.2 X-ray Luminosity function estimation: 
XMM-XXL field 


In the case of the XMM-XXL field, only spectroscopically 
confirmed sources in the redshift interval 2 = 3 — 5 are 
used. For that sample there is no need to apply equa¬ 
tion to determine the corresponding AGN space den¬ 
sity. The limitation of that sample however, is that it is 
both X-ray and optical flux limited because of the mag¬ 
nitude limit of the spectroscopic follow-up observations. 
Both cuts need to be accounted for to infer the X-ray lu¬ 
minosity function via equation We do that by exploit¬ 
ing the fact that the 2 > 3 XMM-XXL sample consists 
of powerful [Lx (2— lOkeV) > 10'*'‘ergs“^] broad-line QSOs. 
We use the well-established correlation between monochro¬ 
matic X-ray [L,/(2keV)] and UV [Li^(2500A)] luminosities 
of broad-line QSOs (e.g. [Steffen et al.|[2006{ [Just et ‘ah] 
2007 Lusso et al.||2010| ) to link the observed optical mag¬ 
nitudes and X-ray fluxes of the sample and account for the 
selection effects. Figure [^demonstrates the correlation be¬ 
tween Lu (2 keV) vs L^ (2500A) using X-ray selected broad¬ 
line QSOs from the XMM-XXL field in the redshift interval 
2 = 1 — 5. The low redshift cut is to avoid X-ray AGN 
with relatively low luminosities, for which the UV/optical 
continuum shows non-negligible contribution from the host 
galaxy. We adopt the [Lusso et al.|(2010| bisector best-fit 
logLj.(2keV) = 0.760 logL„(2500 A)T3.508. At a given 
monochromatic optical luminosity we assume that the data 
points scatter around the above relation following a Gaus¬ 
sian distribution with standard deviation a. Figure shows 
that this is a reasonable assumption. From that figure we 
estimate o = 0.4. Broad Absorption Line (BAL) QSOs rep¬ 
resent up to about 26% of optically selected samples (Hewett 


fc Foltz|2003 Reichard et al.|2003 Gibson et al.|20d9 l, and 

are known to be X-ray faint either because of absorption 
or intrinsic X-ray weakness (e.g. [Gallagher et al.||2006 Luo 
et al.|[2014 |. Such sources do not appear to skew the distri¬ 
butions plotted in Figures [^ and [^ This is likely related to 
the bright X-ray flux limit of the XMM-XXL survey, which 
selects against X-ray faint Under assumptions on the optical 
Spectral Energy Distribution of QSOs, Figure[^can be used 
to estimate the SDSS r-band optical magnitude distribution 
of X-ray sources of given Lx and 2 and then determine the 
fraction of this distribution that is brighter than the spec¬ 
troscopic magnitude limit of the survey, r^ut = 22.5. In this 
case the expected number of detected sources with the sur¬ 
veyed area. A, in equation can be rewritten 


(4) 


A=y dlogLx ^d2 A(Lx, 2) B{Lx,z, \ r) 

X (t){Lx,z I 9) r]{r), 

where B{Lx,z \ r) is the SDSS r-band magnitude distri- 
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bution of a source with Lx and 2 . The efficiency factor 
ri{r) is the success rate of measuring secure redshifts for 
X-ray sources as a function of the SDSS r-band magni¬ 
tude. It accounts for the fact that not all X-ray sources 
with secure counterparts have successful spectroscopic red- 
shift measurements. Collisions between SDSS fibers or the 
finite number of science fibers on the SDSS spectroscopic 
plates mean that not all candidate sources for follow-up 
spectroscopy can be assigned a fiber. Moreover, for the 
sources that are observed the rate of secure redshift mea¬ 
surement depends on their optical brightness. At fainter 
magnitudes the signal-to-noise ratio of the optical spectra 
decreases and therefore the ability to estimate redshifts is 
affected. The probability of a source being assigned a fiber 
is random, while the redshift success rate depends, at least 
to the first approximation, on optical magnitude. The fac¬ 
tor ri{r) is the number of spectroscopically confirmed X-ray 
sources in the magnitude interval r ± Ar divided by the 
total number of X-ray sources that are potential targets 
for follow-up spectroscopy in the same magnitude range. 
The parent X-ray sample of potential targets is selected to 
have /x(0.5 — 10keV,r = 1.4) > ergs“^ cm“^ and 

15 < r < 22.5 mag (see Section [2.2[ ). The r-band magnitude 
dependence of r](r) is shown in Figure]^ For magnitudes in 
the range r = 17.0 — 21.5 mag the efficiency factor r]{r) is 
nearly constant and larger than 80%. This fraction drops to 
about 50% at r = 22.5, the limiting magnitude for follow¬ 
up spectroscopy with the Sloan telescope, and is zero for 
r > 22.5 mag. 

In equation for the calculation of the optical k- 
corrections we use the simulated QSO SEDs of |McGreer| 


et al. (20131. They are generated assuming a double power- 


law continuum in the rest-frame UV/optical part of the 
SED with a break-point at 1100 A. The short- and long- 
wavelength slopes are drawn from normal distributions with 
means of —1.7 and —0.5, respectively. For both slopes the 
scatter is fixed to 0.3. Emission lines with luminosity depen¬ 
dent equivalent widths as well as Lya forest absorption are 
also added to the simulated SEDs as described in IMcGreerl 
et al. ] poT^ . A total of 60 000 model SEDs are generated 


in the redshift interval 2 = 3 — 5 and for AGN luminosities 
logL,y(2500A) ~ 29 — 32ergs“^. These are used to calcu¬ 
late the expected distribution of the observed r-band optical 
magnitudes for a given redshift and AGN luminosity, and de¬ 
termine the term B{Lx,z, \ r) in equation equation]^ For 
the calculation of X-ray k-corrections we assume a power- 
law X-ray spectrum with F = 1.9. 


3.3 X-ray luminosity function models 

The Bayesian framework outlined above explicitly requires 
a model with a set of free parameters that are constrained 
by the observations. We consider both non-parametric and 
parametric models for the X-ray luminosity function in the 
redshift range 2 = 3 — 5. The non-parametric model simply 
assumes that the space density of AGN is constant within a 
given luminosity and redshift interval. In this case the free 
model parameters to be estimated are the AGN space den¬ 
sities in each luminosity and redshift bin. This is equivalent 


to the widely used 1 / Vma>j^pproach for the determination 
of binned luminosity functions. The advantage of the non- 
parametric approach is that it allows investigation of the 
form and amplitude of the X-ray luminosity function evolu¬ 
tion in a model-independent way. 

We also consider four parametric models for the X- 
ray luminosity function and its redshift evolution, which 
have been extensively used in the literature. These are the 
Pure Luminosity Evolution (PLE), Pure Density Evolution 
(PDE), Luminosity Dependent Density Evolution (LDDE) 
and Luminosity And Density Evolution (LADE) models. 
The parametrisation of each model follows [Vito et ah (20141. 
The X-ray luminosity function in the redshift range 2 = 3 — 5 
is defined as the space density of AGN per logarithmic lu¬ 
minosity bin and is described by a double power-law of the 
form 


where K is the normalization, 71 and 72 are the faint and 
bright-end power-law slopes, respectively, and L* is the 
break luminosity. The PLE model assumes that L* is a func¬ 
tion of redshift and evolves according to the relation 


LFz) = LFzo) X , ( 6 ) 

where we fix 20 = 3 and the exponent p parametrises how 
fast the break luminosity evolves with redshift. In the case 
of PDE it is assumed that only the normalisation of the 
luminosity function evolves with the redshift according to 
the relation 


K(.| = K(») X (i^) 


(7) 


where 20 = 3 and the exponent q parametrises the speed 
of the normalisation factor evolution. The LADE model 
adopted here is similar to the Independent Luminosity and 


Density Evolution (ILDE) model described by Yencho et al. 

(20091. We use equation to parametrise the redshift evo¬ 
lution of L* and also add a normalisation evolution term of 
the form 


K{z) = K{zo) X 


( 8 ) 


Finally we consider the LDDE parametrisation of [Hasinger 
et al. (20051, where the normalisation factor of equation [5 


changes with redshift as 


K{z) = Kizo) X 


1 -|- 2 \ 9 +^ -44) 

1 + Zo / 


(9) 


^ In the case of a step function (i.e. non-parametric) model for 
<:l>(Lx,z) and vanishing uncertainties for the redshifts and lumi¬ 
nosities of individual sources, it is straightforward to show that 
the maximum likelihood value of (j){Lx’,z) in equation 1 reduces 
to the [Page &: Carrera] l [2000| | binned luminosity function estima¬ 
tor. 
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Figure 8. AGN space density as a function of 2-10 keV X-ray luminosity in the redshift intervals z = 3 — A and 2 : = 4 — 5. In both panels 
the datapoints are the non-parametric binned X-ray luminosity function. The errors correspond to the 68th percentiles of the Probability 
Distribution Function. The black filled circles are the estimates for the combined XMM-XXL and Chandra deep field surveys. In the 
2 = 3 — 4 panel (left) we also plot the constraints obtained separately from the XMM-XXL (red squares) and the Chandra deep survey 
fields (blue triangles). This is to demonstrate that in overlapping luminosity bins the XMM-XXL and Chandra data yield consistent 
results. The shaded regions are the 68% confidence intervals of the LDDE and PDE parametric models described in Section |3.3| The 
gray hatched region is for the LDDE and the pink region is for the PDE. The parametric models are estimated at the middle redshift of 
the redshift intervals, i.e. 2 = 3.5 and 2 = 4.5. 


where as previously 20 = 3 and the rate of the density evolu¬ 
tion also depends on the X-ray luminosity via the parameter 

/3. 

Recall that in addition to the model parameters that 
describe the X-ray luminosity function in the redshift range 
2 = 3 — 5 the likelihood function in equation]^ also includes 
terms that correspond to the total AGN space density out¬ 
side the redshift range of interest, 2 < 3 or 2 > 5 (see equa¬ 
tion]^. A non-parametric model is adopted to describe that 
term. This to minimise the impact of parametric model as¬ 
sumptions to the results and conclusions. We use 3 redshift 
bins in the range 2 = 0 — 3 with size A 2 = 1. Each redshift 
bin is split into 5 logarithmic luminosity bins in the range 
logLx(2 — lOkeV) =41 — 46 (units of ergs“^) with width 
A log Lx = 1 dex. One additional term is used to model the 
luminosity function in the range 2 = 5 — 6. A constant AGN 
space density is assumed within the above luminosity and 
redshift intervals. We therefore use a total of 16 nuisance 
parameters to model the AGN X-ray luminosity function at 
2 < 3 or 2 > 5. 

The MultiNest multimodal nested sampling algorithm 
I Feroz & Hobson 2008 Feroz et al. 20091 is used for both 
Bayesian parameter estimation and the calculation of the 
Bayesian evidence, Z, of each model, i.e. the integral of the 
model likelihood over the parameter space allowed by the 
priors. The Bayesian evidence is used for model comparison, 
i.e. to select among different models the one that better 
describes the data. 


3.4 Potential obscuration effects 

The analysis presented in this paper uses the full-band se¬ 
lected sample to determine the X-ray luminosity function 
of AGN at 2 > 3. This is because of the higher sensitivity 


of the 0.5-10 keV band compared to e.g. the 2-10 keV band, 
which translates to a larger number of sources. A potential 
issue however, is that the analysis described above ignores 
the impact of obscuration on the observed AGN flux. Fig¬ 
ure shows that for AGN at 2 > 3 moderate line-of-sight 
columns, log Ah = 23 (cm“^), suppress the observed 0.5- 
lOkeV flux of AGN by less than about 25%. Higher column 
densities however, have a larger impact on the 0.5-10 keV 
flux and therefore result in incompleteness in the sample 
that is not accounted for in our analysis. Nevertheless, we 
are primarily interested in the differential evolution of AGN 
in the redshift intervals 2 = 3 — 4 and 2 = 4 — 5. Figure]^ 
shows that for hxed obscuration and intrinsic AGN luminos¬ 
ity there is little difference in the level of flux suppression 
between redshifts 2 = 3 and 2 = 5. Incompleteness related 
to obscuration is therefore expected to be similar across the 
redshift range 2 = 3 — 5. This allows direct comparison of the 
inferred AGN space densities in the redshift bins 2 = 3 — 4 
and 2 = 4 — 5, under the assumption that the distribution of 
AGN in obscuration does not change dramatically between 
these redshift intervals. Recent studies on the obscuration 
distribution of AGN (Ueda et al.|2014 Buchner et al.|2015 


|Aird et al.||2015p show that the obscured AGN fraction in¬ 
creases with redshift, at least to 2 « 3. In [Buchner et ah] 
(20151 for example, the obscured fraction in Gompton thin 
AGN is about 50% at 2 < 1 and increases to 70% at 2 « 3. 
At higher redshifts however, current constraints on the evo¬ 
lution of the obscured AGN fraction are still limited by small 
number statistics. Nevertheless, there are suggestions that 
the obscured AGN fraction remains roughly constant with 
redshift at 2 > 3 | Buchner et al.|2015 I. 


Obscuration related effects potentially have a larger 
impact on the analysis of the XMM-XXL data. For these 
sources it is explicitly assumed that their optical/UV con- 
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Figure 9. X-ray luminosity dependence of the ratio of the AGN 
space densities in the redshift intervals 2 : = 3 — 4 and z = 4 — 5. 
The data points are for the non-parametric binned X-ray lumi¬ 
nosity function. A double power-law is also fit independently to 
the data in the redshift intervals z = 3 — 4 and z = 4 — 5. The ratio 
of these parametric fits is shown by the light/dark grey shaded re¬ 
gions. The widths of the dark-grey and light-grey shaded regions 
at a given luminosity corresponds to the 68% and 95% confi¬ 
dence intervals. The predictions of the Gilli et al. (2007) LDDE 
parametrisation of the luminosity function with (dashed line) and 
without (dot-dashed line) an exponential cutoff at high redshift 
is also plotted for comparison. 

tinua are not obscured or reddened. This assumption allows 
us to place them on the Li, (2 keV) — Li, (2500 A) correlation 
for type-I AGN and correct for the spectroscopic magni¬ 
tude cut as explained above. This poses a problem when 
comparing the Chandra with the XMM-XXL constraints 
on the X-ray luminosity function. Nevertheless, the shal¬ 
low XMM-XXL survey is only sensitive to powerful QSOs 
[Lx (2 —10 keV) > 10'*^ erg s”'^] at 2 > 3. At lower redshift at 
least, it is well established that the obscured AGN fraction 
is rapidly decreasing with increasing accretion luminosity 


(e.g. Ueda et aL||2003 Akylas et al. 

2006 

Ueda et aL||2014 

Buchner et al. 2014| Merloni et al. 

2014 

1 . Obscuration in- 


completeness corrections are therefore likely to be small for 
luminosities Lx(2 — lOkeV) > 10"^^ergs“^. We explore this 
assumption further in the next section by directly compar¬ 
ing the Chandra and XMM-XXL constraints in overlapping 
luminosity bins. 


Table 2. Power-law fits to the data in the redshift intervals z = 
3 — 4 and 2 = 4 — 5 


redshift 

log A 

logL* 

71 

72 

interval 

(Mpc“5) 

(erg/s) 



( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

2 = 3-4 

2 = 4-5 

0-iO_o.i7 

D.OZ_q g2 

44 41 +013 

45 00 

^O.UU_o_53 

0.l9t°;l®3 

0-68tJ!:"38 

9 95+0-29 

2 16+0-98 

^•lO_o.47 


Listed are the best-fit parameters for the simple power-law fits 
to the data in the redshift intervals 2 = 3 — 4, z = 4—5. The listed 
values are the median of the probability distribution function of 
each parameter. The errors correspond to the 16th and 84th per¬ 
centiles around the median. The columns are: (1) redshift inter¬ 
val (2) X-ray luminosity function normalisation (see equation]^, 
(3) break luminosity of the X-ray luminosity function (see equa¬ 
tion]^, (4) faint-end slope, (5) bright-end slope. 

lOkeV) = 42.5 (ergs“^, but do not probe sufficient vol¬ 
ume to detect luminous and rare sources with luminosities 
> 10"^® ergs“^. The XMM-XXL survey can constrain the 
bright-end of the luminosity function but is not sensitive 
enough to detect AGN below Lx (2 — 10 keV) < 10^'* erg s“^. 
Nevertheless, in the interval log Lx (2 — 10 keV) =44 — 45 
(ergs“^) both the Chandra and XMM data provide mean¬ 
ingful constraints to the AGN space density. These indepen¬ 
dent non-parametric estimates are in good agreement within 
the 68 % errors plotted in Figure This suggests that ob¬ 
scuration effects do not have a strong impact on the AGN 
space densities estimated from the XMM-XXL data. In Fig¬ 
ure l^the 2 = 4 — 5 panel does not plot separately the space 
density estimates from the Chandra and then XMM-XXL 
survey fields. This is because at this redshift interval there 
is virtually no overlap in the luminosities probed by the two 
datasets. The XMM-XXL constrains only the brightest lu¬ 
minosity bin, logI/x(2 — lOkeV) = 45 — 46 (units ergs“^). 
The Chandra fields include only sources that are fainter than 
Lx(2-10keV) = 10'‘®ergs-i . 

A striking result illustrated in Figurej^is the strong evo¬ 
lution of the AGN population between redshifts 2 = 3 and 
2 = 5. This is further demonstrated in Figure]^ which plots 
as a function of luminosity the ratio of the non-parametric 
AGN space densities in the redshift intervals 2 = 3 — 4 and 

2 = 4—5. At luminosities I/x(2 —lOkeV) < 10'*® ergs“^ there 
is a factor of 5 decrease in the AGN space density from 2 = 

3 — 4 to 2 = 4 — 5. There is also evidence for luminosity de¬ 
pendent evolution. AGN with Lx (2 — 10 keV) < 10*® ergs“* 
in Figure [^appear to evolve faster than intrinsically brighter 
sources. We quantify the statistical significance of the evi¬ 
dence for luminosity dependent evolution using the quantity 


4 RESULTS 

4.1 The non-parametric X-ray luminosity 
function determination 

Figure]^ plots the non-parametric estimate of the X-ray lu¬ 
minosity function in two redshift intervals, 2 = 3 — 4 and 
2 = 4 — 5. For the low redshift sub-sample, 2 = 3 —4, we also 
plot separately the AGN space density determined indepen¬ 
dently from the Chandra deep fields and the XMM-XXL sur¬ 
vey. As expected these datasets probe different luminosity 
ranges. The Chandra fields provide constraints to log Lx (2— 


7^ =log 


log 


0(Lx = 10*® - 10*®,2 = 4-5) 
0(Lx = 10« - 1046, 2 = 3-4)' 

0(Lx = 10*® - 10*®,2 = 4-5) 
4>{Lx = 10 « - 1045 , 2 = 3 - 4)’ 


( 10 ) 


where (f>{Lx , 2 ) is the binned (non-parametric) luminos¬ 
ity function. TZ is the logarithmic difference between the 
high luminosity {Lx > 10 *® ergs”®) data point of Figure]^ 
and the sum of the two moderate luminosity data points 
{Lx ~ 10*® — 10*® ergs”®) of the same plot. We bin to¬ 
gether moderate luminosity AGN to simplify the problem 
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and also improve the statistical reliability of the results. We 
use the full probability density distribution of TZ and find 
that at the 90% probability 7?. > 0, i.e. there is differential 
evolution between moderate {Lx = 10 ^^ — 10 "^® ergs“^) and 
powerful {Lx > 10'*® ergs”*) X-ray AGN from 2 = 3 — 4 to 
2 = 4 — 5. We further investigate this using a simple para¬ 
metric approach. A double power-law function (equation]^ 
is fit independently to the data in the redshift bins 2 = 3 — 4 
and 2 = 4 — 5. In this exercise evolutionary effects within 
each of the two redshift intervals are ignored. The best-ht 
parameters are presented in Table The ratio between the 
two double power-laws in the redshift bins 2 = 3 — 4 and 
2 = 4 — 5 is plotted with the shaded region in Figure 
The apparent increase of this ratio toward faint luminosi¬ 
ties is because the faint-end slope of the 2 = 4 — 5 sam¬ 
ple is poorly constrained and on the average steeper than 
that of the 2 = 3 — 4 sample. Nevertheless, in the interval 
Lx(2 — lOkeV) = 10*® — 10*® ergs”*, where constraints on 
the AGN space density are available for both the 2 = 3 — 4 
and the 2 = 4—5 sub-samples, the ratio between the two 
power-laws is consistent with luminosity-dependent evolu¬ 
tion. Larger samples are needed however, particularly in the 
redshift interval 2 = 4 — 5, to reduce the uncertainties in 
Figure and further explore the evidence for luminosity de¬ 
pendent evolution. 

The amplitude of the AGN X-ray luminosity function 
evolution between the redshift intervals 2 = 3 — 4 and 
2 = 4 — 5 in Figure is further compared with the predic¬ 
tions of luminosity function parametrisations that include 
an exponential decline at 2>3 (e.g. Gilli et al.||2007l. Such 
models are motivated by the rapid evolution ot optical QSO 


space density at hig h redshift (e.g. Schmidt et aJ. 1995 
Richards et al.|2006l. We use the Gilli et al. (2007)[^jLDDE 


parametrisation of the X-ray luminosity function with and 
without an additional exponential cutoff to predict the space 
density of AGN at 2 = 3.5 and 2 = 4.5. The ratio be¬ 
tween these predictions is plotted in Figure The ampli¬ 
tude of the AGN evolution inferred in this paper for lu¬ 
minosities Lx(2 — lOkeV) < 10*® ergs”* is consistent with 
the Gilli et al. (2007) LDDE parametrisation that includes 
an exponential cutoff. AGN with luminosities in the range 
Lx (2 — 10 keV) = 10*® — 10*® erg s”* lie in between the Gilli 
et al. (2007) model predictions with and without an expo¬ 
nential cutoff. This is suggestive of milder evolution, albeit 
at the « 90% confidence level. 


4.2 Parametric X-ray luminosity function 
determination 

Next we use Bayesian model comparison to assess which of 
the evolutionary models outlined in Section |3.3| provides a 
better description of the observations in the redshift interval 
2 = 3 — 5. The PLE, PDE, LADE and LDDE parametric 
models of Section [3.31 are ht to the combined Chandra and 
XMM-XXL dataset. Table presents the best-fit parame¬ 
ters for each parametric model. The Bayes factor (ratio be¬ 
tween evidences) of each model relative to the model with 
the highest evidence (PDE) is also shown in that table. 

The model with the highest evidence in Table is the 


^ http://www.bo.astro.it/-gilli/counts.html 


PDE, with the second best being the LDDE. The Bayes fac¬ 
tor of the two models is A logjQ Z = 0.25. Based on the 
Jeffreys interpretation of the Bayes factor (Jeffreys 19611 
this difference suggests that both models describe equally 
well the evolution in the redshift interval 2 = 3 — 5 of the X- 
ray selected AGN sample presented in this paper. The bulk 
of the AGN in the present sample have X-ray luminosities 
Lx(2 — lOkeV) < 10*® ergs”*. Figure § shows that such 
sources experience similar reduction in their space density 
between 2 = 3 — 4 and 2 = 4 — 5, i.e. consistent with pure 
density evolution. More luminous AGN [Lx (2 — lOkeV) > 
10 *® erg s”*], which appear to experience milder evolution in 
Figure[^ represent only a small fraction of the present sam¬ 
ple. This combined with the fact that the PDE is a simpler 
model compared to the LDDE (5 vs 6 free parameters) re¬ 
sults in similar Bayesian evidences for these two parametric 
models. 

Table also shows that the LADE parametrisation 
adopted in this work performs worse than the LDDE. The 
Bayes factor of the two models is Alog^gZ = 1.08. This 
difference is strong evidence in favor of the LDDE |Jeffreys| 
| T96T| ). 

The one model that performs significantly worse than 
the rest is the PLE. The Bayes factor of that model relative 
to the one with the highest evidence (PDE) is A logj^Q Z = 
4.02. Based on the Jeffreys interpretation of the Bayes fac¬ 
tor | |Jeffreys|1961[ ) this is decisive evidence against the PLE 
model. 

Although the PDE model is favoured by our analysis 
for the evolution of AGN in the redshift interval 2 = 3 — 5 in 
the next sections we use the LDDE model to compare with 
previous studies and to determine the contribution of X-ray 
AGN to the ionisation of the Universe. This is because of the 
small difference in the evidences of the PDE vs the LDDE 
and the long literature on the LDDE evolutionary model. 
Nevertheless, the results and conclusions are not sensitive 
to the particular choice of AGN evolutionary model. 


4.3 Comparison with previous studies 


A number of studies on the X-ray luminosity function of 
AGN have appeared in the literature recently. These include 
works focusing on the space density of X-ray AGN at high 


redshift s 2 : > 3 ( 

Civano et al. 12010 

Kalfountzou et al. 

2014 

Vito et al.|2014 

1 and results on the global X-ray luminosity 


function evolution across redshi ft (|Ueda et al.||2014[ |Buch-| 
ner et al.|2015 Aird et al.|2015 Miyaji et al.|2015 |. In this 

section we compare our X-ray luminosity function with pre¬ 
vious studies using X-ray samples selected in the 0.5-2 or 0.5- 
lOkeV bands, i.e. similar to the one presented in this paper. 


Vito et al. (20141 compiled one of the largest samples to date 


in the redshift interval 2 = 3 — 5. They select sources in the 
0.5-2 keV band and also present non-parametric (1/Vmax) 
estimates of the X-ray luminosity function. These aspects of 
the Vito et al. analysis methodology are similar to ours. We 


also use recent results of Aird et al. (20151 as a represen¬ 


tation of parametric approaches to fit the X-ray luminosity 
function of AGN to the full redshift interval accessible to 
current X-ray selected samples, z ~ 0 — 5. 

Vito et al. (2014) use a sample similar to ours in size 
to determine the X-ray luminosity function in the range 
3 < 2 < 5. They use both photometric and spectroscopic 
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redshifts (total of 141) and apply corrections to account for 
sources without photometric redshift determinations (total 
of 65). We compare our best-fit LDDE model with their re¬ 
sults in Figure |10| The redshift intervals of each panel are 


the same as in Figure 7 of Vito et al. (2014|. We plot their 
non-parametric 1 /Vmax estimates that include their redshift 
incompleteness corrections. Our LDDE parametric models 
in Figure is estimated at the middle of each redshift in¬ 
terval. We find that our X-ray luminosity function determi¬ 
nation is systematically lower than the Vito et al. (2014) 
data points. This discrepancy (> 3cr significance in e.g. the 
redshift interval z = 3.47 —3.90 of Figure [To| ) may be related 
to the fact that Vito et al. (2014) account for X-ray obscu¬ 
ration in the determination of X-ray luminosities and the 
calculation of the Vmax of individual sources. A total of 36 
of their 141 sources have column densities Nh 10 ^® cm“^, 
which translates to suppression of their fluxes compared to 
the unobscured (Nh = 0cm“^) case by >20% (see Fig¬ 
ure [^. Such corrections are ignored in the analysis pre¬ 
sented here. Alternatively the discrepancy may be related 
to how sources with photometric redshifts or sources with¬ 


out any redshift information are treated. Vito et al.| (20141 


use only the best-fit photometric redshift solution without 
taking into account the corresponding uncertainties. It is 
further assumed that the 65 sources in their sample without 
photometric redshifts all lie in the redshift interval 2 = 3 — 5 
and that they follow the redshift distribution of the X-ray 
sources with redshift determinations (photometric or spec¬ 
troscopic) in the range 2 = 3 — 5. The amplitude of this 
correction is « -1-0.5 and -1-0.25 dex increase at luminosities 
Lx(2 — lOkeV) « lO'^^ and 10"^® ergs“^ respectively. In that 
respect [Vito et ^ (20141 determine maximal X-ray lumi¬ 
nosity functions. It is likely that some of the X-ray sources 
without redshift information lie outside the redshift range 
2 = 3 — 5. Heavily obscured AGN or AGN hosted by early 
type hosts at moderate redshifts, 2 « 1 — 3, may also have 


red SEDs, similar to those of high redshift sources ( 

Koeke- 

moer et al.|2004 

[Schaerer et al.|2007||Rodighiero et 

al. [20071 

Del Moro et al. 

2009 



We further investigate this issue by adopting a method¬ 


ology similar to that of Vito et al. (20141 to determine the 


X-ray luminosity function. We use the best-fit photometric 
redshifts only, i.e. ignoring the photometric redshift proba¬ 
bility distribution functions. Sources without optical iden¬ 
tifications in the sample are assigned random redshifts in 
the range 2 = 3 — 5 based on the redshift distribution of 
sources with spectroscopic or photometric redshift measure¬ 
ments. The results are presented in Figure It shows that 
an approach that ignores photometric redshift errors results 
in an overestimation of the AGN space density. This is not 
only because of sources without optical identifications and 
therefore without photometric redshift (total of 107 in our 
sample) estimates being forced to lie in the range 2 = 3 — 5. 
The overestimation is mainly the result of ignoring photo¬ 
metric redshift uncertainties. In the likelihood equation[^the 
probability of a source having redshift 2 and luminosity Lx 
is weighted by the luminosity function. Sources with broad 
(see Fig. redshift probability distribution functions are 
more likely to lie at low redshift and moderate luminosities 
simply because the space density of AGN is higher there. 
Ignoring this weighting in equation by e.g. fixing photo¬ 


metric redshifts to a single (best-fit) value overestimates the 
luminosity function at high redshifts. 

Also plotted in Figur e|10|are the flexible double power- 
law parametric model of Aird et al. (20151 for their 0.5-2 


selected sample before applying corrections for obscuration 
(i.e. directly comparable to our analysis). The binned X-ray 
luminosity function estimates of Aird et al. (20151 are also 


shown. These data points are estimated using the Nobs/Nmdi 


method developed by Miyaji et al. (20011 and are therefore 


not independent of the underlying parametric model plotted 
in Figure Nevertheless, these binned estimates provide 
a measure of the associated uncertainties and can be used 
to identify redshift and luminosity intervals where the para¬ 
metric model provides a poorer fit to the data. The |Aird| 
et al. (20151 X-ray luminosity function estimated from their 


0.5-2 keV selected sample is in fair agreement with our re¬ 
sults. This may be not surprising given the overlap between 
the two datasets and the similar approaches. Figureplots 
the X-ray luminosity density evolution of AGN using our 


LDDE model and the Aird et al. (20151 total X-ray lumi¬ 


nosity function that includes both obscured (Compton thin 
and thick) and unobscured sources. The X-ray luminosity 
function estimated in this work accounts for about 40% of 
the total X-ray luminosity density determined by|Aird et al.| 


Finally, we also compare in Figure[^the X-ray luminos¬ 
ity function with high redshift determinations of the optical 
QSO luminosity function. The conversion from UV/optical 
to X-rays ultimately depends on the scatter in the relations 
between bolometric and X-ray or UV luminosities. There 
are suggestions that the bolometric-to-X-ray luminosity ra¬ 
tio has a larger scatter than the bolometric-to-UV luminos¬ 
ity ratio ( [Hopkins et al.|[2007[ ). We therefore convert from 
UV/optical to X-rays by convolving the best-fit parametric 
model of UV/optical QSOs from the relevant publications 
with the Ly(2keV) — Li,(2500A) relation of 


Lusso et al. 


(20101 assuming a scatter of 0.4 dex (e.g. Figure [^. This 


calculation also requires assumptions on the shape of the 
UV SED of QSOs at wavelengths A > 1450 A. We assume 
a power-law of the form L{v) oc v~^'^ (Vanden Berk et al. 


2001 Teller et al. 20021. A steeper slope, L{i/) oc 
( [Lusso et al.|2015 l, translates to an X-ray luminosity func¬ 
tion that is offset by 5 log Lx ~ —0.03 dex compared to 
L{v) oc v~^'^. At the 2 = 3.1 and 2 = 4.1 panels of Figure [l0[ 
we overplot the Masters et al. (20121 best-fit double power- 


law models at 2 = 3.2 and 2 = 4. At the 2 = 3.35 panel of 
Figure 10 we show the [Ross et al.[ ( |2013| LEDE (L uminos- 
ity Evolution and Density Evolution) model fit to the BOSS 
Stripe82 QSO data ( [Palanque-Delabrouille et al.|20lT| . We 
choose not to extrapolate this model past redshift 2 = 3.5, 
the limiting redshift of the BOSS Stripe82 QSO sample. At 
the redshift interval 4.3 < 2 < 5.1 (mean redshift 2 = 4.7) 
we transform to X-rays the optical QSO luminosity function 
determined by McGreer et ah] (20131. The optical luminosity 
functions are plotted by shaded regions in Figure [10[ The 
shape and normali sati on of the UV/optical QSO luminosity 
functions in Figure 10 at luminosities > 10'*® ergs”* are sen¬ 


sitive to the scatter of the L„ (2 keV) — (2500 A) relation 

used in the convolution. 

A striking result from Figure is the steep faint-end 


slope of the Masters et al. (20121 best-ht double power-law 
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models at 2 = 3.2 and 2 = 4 compared to the faint-end slope 
of the X-ray luminosity function. The ratio between the 
UV/optical and the X-ray luminosity functions is a proxy of 
the type-I fraction among AGN, 7uype-i, and can be used to 
explore the luminosity dependence of this quantity at high 
redshift. In this section we define type-I AGN as sources 
with blue UV/optical continua, typical of broad-line QSOs. 
We use the Masters et al. (2012) luminosity function and 
the LDDE parametrisation for the X-ray luminosity func¬ 
tion to determine 7uype-i and plot the results in the case 
of 2 = 3.2 in Figure |13| The luminosity dependence of the 
Jlype-i at 2 = 4 is very similar and is not shown. We find 
evidence that the type-I AGN fraction at 2 > 3 is a non¬ 
monotonic function of luminosity. There is a minimum at 
Lx(2 — lOkeV) « lO'^^ergs”^, followed by a steep increase 
toward fainter luminosities. 

The behaviour of Jlype-i in Figure is opposite to 
studies that find a drop in the type-I or X-ray unobscured 
(Nh < 10^^ cm“^) AGN fraction with decreasing luminosity 
or equivalently that type-II or obscured (Nh > 10^^cm~^) 


AGN dominate at faint luminosities (Ueda et al. 12003 


las et alT]|2006| |Ueda et al.||2014{ [Merloni et al.||2014 


Aky 


et al.|20i5 I. Studies on the obscuration distribution of AGN 


in the local Universe support a picture where the obscured 
AGN fraction increases with decreasing luminosity. Never¬ 
theless, they also find evidence for a turnover (drop) of the 
obscured AGN fraction at v ery faint X-ray luminosities, be- 
low about lO^^ergs”^ (e.g Burlon et al. 2011 Brightman 


fc Nandra|2011b I. Recently, Buchner et al. (20151 extended 


these results to higher redshift and found evidence that the 
obscured AGN fraction peaks at a redshift-dependent lumi¬ 
nosity and then drops at both brighter and fainter luminosi¬ 
ties. Figure [T3| overplots the obscured AGN fraction derived 
by Buchner et al. (2015) in the redshift interval 2 = 2.7 — 4. 
For this comparison we assume that Luype-i = 1 — L^obscured. 
This is a simplistic assumption because the definition of 
Type-I QSOs in the case of UV/optical samples and un¬ 
obscured/obscured AGN in X-ray samples like in Buchner 
et al. (2015) is different. Nevertheless, to the first order the 
quantity 1 — J4)bscured should be at least loosely related to 
Juype-i. In Figure fT^ there is qualitative agreement between 
the Buchner et al. (2015) 1 — J^obscured parameter and our 
definition of Juype-i. 

At bright luminosities the quantity Juype-i in Fig¬ 
ure [^increases with increasing Lx - This in agreement with 


previous studies based on X-ray (|Ueda et al. 

|2014| e.g.). 

optical (e.g. Simpson 20051 or infrared (e.g. 

Assef et al. 


2013) data. At the brightest luminosities probed by our 
data. Lx = 10"^® — 10^®ergs“^, the type-I fraction is about 
75 ±25%. This number is relevant to the population of pow¬ 
erful QSOs (bolometric Lboi ~ 10"^^ergs“^) with reddened 
UV/optical continua [extinction F(B — V) « 5] identified in 
recent wide-area infrared surveys ( |Stern et al.|2014 |. These 
reddened/obscured sources corresp ond to X-ray luminos ities 
Lx(2 — 10keV) ~ 5 x 10"^® ergs“^ ( Marconi et al.|2004 ) and 
are suggested to represent up to 50% of the QSO popula¬ 
tion at these luminosities (Assef et al. 2014). We caution 
however, that there are uncertainties on the inferred bolo¬ 
metric luminosities of these sources and hence, the obscured 
AGN fraction, depending on the assumed geometry, physical 


scale and covering fraction of the obscuring material (Assef 
et al.|[2014|. Additionally, the X-ray properties of these in- 



Figure 12. X-ray luminosity density as a function of redshift. 
The black dashed curve corresponds to the total X-ray luminos¬ 
ity density for both obscured (Compton thin and Compton thick) 
and unobscured AGN deterined by Aird et al. (2015). The shaded 
blue region is the X-ray luminosity density estimated using our 
LDDE parametrisation for the evolution of the AGN population 
at 2 > 3. The inset plot shows the ratio between the blue (solid) 
and black (dashed) curves as a function of redshift. It represents 
the fraction of the luminosity density accounted by our LDDE 
model relative to the total luminosity density (corrected for ob¬ 
scuration effects) estimated by the Aird et al. (2015). 


frared selected AGN are poorly known. There are sugges¬ 
tions that they represent heavily obscured, even possibly 
Gompton thick, systems (Stern et al. 2014), in which case 
they are expected to be underrepresented in our sample. 


4.4 Contribution of AGN to ionisation of the 
Universe 

The X-ray luminosity function can be used to set limits 
on the contribution of the AGN to the radiation field that 
keeps the Universe ionised at redshift 2 > 3. The advan¬ 
tage of X-ray selection is that it provides a better handle 
on the faint-end of the AGN luminosity function compared 
to UV/optically selected samples. The downside is that as¬ 
sumptions on the escape fraction of UV photons have to 
be made to convert AGN space densities to ionising photon 
densities. At the very least however. X-rays surveys can set 
strict upper limits on the contribution of AGN to the ion¬ 
ising photon field, under the assumption that all photons 
emitted by the central source escape to the Inter-Galactic 
Medium. 

The rate of hydrogen ionising photons is estimated by 
integrating the AGN Spectral Energy Distribution in the 
energy range 1-4 ryd 


r 

J Iryd 


L(v) 

hv 


dv- 


( 11 ) 


where L(v) is the AGN monochro matic luminosity at fre¬ 
quency V. The Lusso et al. (2010) L,^(2keV) — L,^(2500 A) 
relation is used to convert the 2-lOkeV X-ray luminosity to 
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Figure 10. Comparison of our best-fit LDDE model (pink shaded region) with previous estimates of the X-ray luminosity function in 
the range 2 = 3 — 5. The panels correspond to the redshift intervals of Vito et al. (2014) to allow direct comparison with their results. The 
black filled circles are the 1/Vmax binned luminosity function estimates of Vito et al. (2014). The thin blue curves are the flexible double 
power-law parametric model of Aird et al. (2015) for their 0.5-2 keV band selected sample without corrections for obscuration. For the 
Aird et al. (2015) 0.5-2 keV sample the binned X-ray luminosity function estimates are also shown by the red crosses. The UV/optical 
QSO luminosity functions of Masters et al. (2012), McGreer et al. (2013) and Ross et al. (2013) are also plotted by the thick black lines 
at the relevant redshift intervals. 
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Figure 11. Impact of ignoring photometric redshift uncertainties on the inferred AGN space density. The red open circles correspond 
to the non-parametric binned X-ray luminosity function estimated by using the best-fit photometric redshift solution only and ignoring 
photometric redshift uncertainties (see text for details). The black data points correspond to the AGN space density estimated by using 
the full photometric redshift Probability Distribution Function to account for photometric redshift errors. For both set of data-points 
only the Ghandra survey fields are used. 


UV monochromatic luminosity at 2500 A. We then extrapo¬ 
late to the wavelength range 227A (4ryd) and 910A (Iryd) 
assuming a double power law for the AGN SED of the form 


L{i') oc 


,- 0.5 

,.- 1-7 


(1100 A< A < 2500 A) 
A < 1100 A 


( 12 ) 


The spectral slopes in the above relation are from |Vanden| 
Berk et al. (20011 and Teller et al. (20021. At any given 


redshift the comoving density of the hydrogen ionising rate is 
then estimated by integrating the X-ray luminosity function 


JV = 




X,z) fesc{Lx,z)h dlogLx. 


(13) 


The LDDE parametrisation of the X-ray AGN luminos¬ 
ity function is adopted. The integration limits are set to 
Lx(2 — lOkeV) = 10^^ and 10^®ergs“^. The photon es¬ 
caping factor, feac{Lx, z), accounts for the fraction of ob¬ 
scured AGN, which likely depends on both redshift and ac¬ 
cretion luminosity. In these sources the ionising photons are 
absorbed locally and therefore do not contribute to the cos¬ 
mic ionisation radiation field. Our analysis does not con¬ 
strain the distribution of AGN in obscuration. Additionally 
at redshift z > 3 there are still discrepancies among different 
studies on the obscured AGN fraction and its dependence on 
luminosity (e.g. |Ueda et al.|201^ [Buchner et al.|[2015| jAirdj 
et al.|20l5 I. For example, in Figure 13 we find evidence that 
the type-I AGN fraction and by proxy the obscured AGN 
fraction, have a different dependence on luminosity com¬ 
pared to relations established at lower redshift. We choose 
to present our baseline results for the luminosity dependent 


type-I AGN fraction determined by Merloni et al. (20141. In 


that study type-I refers to AGN with either blue UV/optical 
continua and/or broad emission lines. This definition is more 
releavant to the determination of the hydrogen ionising pho¬ 
ton rates of AGN, compared to e.g. the standrard X-ray un¬ 


obscured AGN dehnition, Nh < 10^^ cm Adopting the 
Merloni et al. (20141 relation for fesc{Lx,z) also allows di¬ 


rect comparison with previous X-ray studies that also used 
monotonically increasing obscured AGN fraction with de¬ 
creasing X-ray luminosity to approximate the escape frac¬ 
tion of UV photons in high redshift AGN. Figure plots 
A/” as a function of redshift for our baseline model for the 
photon escaping fraction. In that figure we also place an up¬ 
per limit in A/” by setting ff.sc = 1, i.e. the extreme case that 
all photons escape and contribute to the hydrogen ionising 
radiation. This translates to a net increase of A/” by a factor 


of about 1.4 compared to the Merloni et al. (20141 type-I 
AGN relation for fssc- We further explore how the results 
in Figure change if fesc is approximated by X-ray defi¬ 
nitions of the unobscured (Nh < 10^^ cm“^) AGN fraction. 
We adopt the luminosity dependence of the unobscured X- 


ray AGN fraction derived by Buchner et al. (20151 in the 


redshift interval a = 2.7 — 4 (see Fig. 131. This assumption 
yields photon rate densities that are a factor of about 1.7 
smaller compared to the baseline results that use the [Mer^ 
loni et al. (20141 type-I AGN fraction as proxy of fesc- For 


clarity these results are not plotted in Figure [T4| 

The constraints above should be compared to the min¬ 
imum photon rate density required to keep the Universe 
ionised at a given redshift. This is estimated from equation 
(26) of jMadau et al. (1999 \ after scaling it to our cosmology. 
We also assume that the ionised hydrogen clumping factor 
in that relation, which is a measure of the inhomogeneity of 
the medium, evolves with redshift as 


C = 1 -|- 43 X 2 


(14) 


This relation is based on cosmological simulations (Pawlik 
et al.|[2009| ) and is adopted by Haardt & Madau (20121 to 
synthesise the evolving spectrum of the diffuse UV/X-ray ra¬ 
diation field. Figure [T4| plots the redshift dependence of the 
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logLj;(2-10 keV) (erg/s) 


Figure 13. Type-I AGN fraction as a function of X-ray luminos¬ 
ity. The red shaded region is the ratio between the z = 3.2 lumi¬ 
nosity functions of UV/optical QSOs from Masters et al. (2012) 
and our LDDE parametrisation for X-ray AGN. The black shaded 
region is 1 — .ToUscured i where .Fobscured is the obscured AGN frac¬ 
tion of Buchner et al. (2015) for the redshift interval z = 2.7 — 4. 


photon rate density reqnired to keep the Universe ionised. 


This figure shows that for the baseline model (i.e. Merloni 


et al. 2014 as proxy of /esc) the AGN contribntion to the 


photon rate density required to keep the Univese ionised 
decreases from 70% at z = 4 to about 20% at z = 5. Assum¬ 
ing /esc = 1 the fractions above translate to npper limits 
< 100% at 2 = 4 and < 30% at 2 = 5. These numbers are in 
broad agreement with some previous studies on the role of 
X-ray AGN in the ionisation of the Universe ([Barger et al.| 


2003a Haardt fc Madau|2012 Grissom et al.|2014 |. The de¬ 

creasing ionising photon rate density fractions with increas¬ 
ing redshift is a direct consequence of the strong evolution of 
the X-ray AGN luminosity function between redshifts 2 = 3 
and 2 = 5. 

Finally in Figure the model constraints on the ion¬ 
ising photon rate density at 2 = 3 — 5 from our X-ray 
sample are compared to previous works based on either 
UV/optical QSO surveys (Glikman et al.|2011 Masters et al. 


2012[ McGreer et al. 20131 or UV/X-ray selected samples 
| |Giallongo et ah 20151. For the calculation of the ionis¬ 
ing photon rate densities the luminosity functions estimated 
in these works are integrated between absolute magnitudes 
M(1450A) = —18 and —28 mag by adopting the UV spec¬ 
trum of Equation |12[ Our baseline model assuming the lMer-] 
|loni et ah] ( |2014[ ) escaping fraction is consistent with the 
lower-normalisation UV/optical data points in Figure 14 


5 DISCUSSION 

In this paper we explore the evolution of the AGN X-ray 
luminosity function in the redshift interval 2 = 3 — 5 by 
combining deep Ghandra surveys with the wide-area and 
shallow XMM-XXL sample. This dataset provides sufficient 
depth and volume to determine the space density of AGN 
over 3dex in luminosity, log Lx (2 — lOkeV) Ri 43 — 46, in 


Figure 14. Hydrogen Ionising photon rate density as a func¬ 
tion of redshift. The shaded regions are the constraints from our 
analysis using the LDDE parametrisation for the X-ray lumi¬ 
nosity function and under different assumptions on the escape 
fraction of AGN photons. The grey-shaded region assumes an 
escaping fraction of unity, i.e. ignoring obscuration effects close 
to the supermassive black hole. The pink-shaded region assumes 
the luminosity-dependent Type-1 AGN fraction of Merloni et al. 
(2014). We caution that beyond z = 5 the shaded curves are ex¬ 
trapolations. The data points correspond to results in literature. 
The thick black line in the plot shows the photon rate density 
required to keep the Universe ionised at any given redshift. The 
ratio between the shaded regions and the black line are presented 
in the inset plot. 



logLx(2-10 keV) (erg/s) 


Figure 15. Number of high redshift AGN as a function of 2- 
10 keV X-ray luminosity that the eROSITA 4-year All Sky Survey 
is expected to detect. The LDDE parametrisation of the X-ray lu¬ 
minosity function is used for the predictions. Results for 3 redshift 
intervals are plotted, 2 = 3 — 4 (blue solid), 2 = 4 — 5 (red dashed) 
and 2 = 5 — 6 (black dot-dashed). The predictions for the latter 
redshift bin are extrapolations of the LDDE model. The numbers 
are for 0.5 dex wide luminosity bins. 
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Table 3. Parametric model best-fit parameters and uncertainties 


Model 

(1) 

log K 
(Mpc-3) 

(2) 

logL* 

(erg/s) 

(3) 

7i 

(4) 

72 

(5) 

p 

(6) 

g 

(7) 

0 

(8) 

logio 2: 

(9) 

A lOg^g ^ 

(10) 

PDE 

LDDE 

LADE 

PLE 

-4.840l°;;| 

-4.79l“;14 

-5.05l“;i4 

44.33l°;13 

44.27l°;lJ 

44.42l[;;l« 

o-2ilg:ll 

o.i8l°:l® 

0-19lal3 

Z.UO_Q 20 

9 it::+0.24 

^•-Lo_0.21 

2 10+0-26 
^•-*-^-0.22 

1 7*3+019 

-0.85l°;;^ 

-4.98t°;«° 

-6.59l°;|° 

-7-4611;“ 

2-2511:15 

2-30ll;|° 

9894.95 

9894.70 

9893.62 

9890.93 

0.0 

-0.25 

-1.33 

-4.02 


Listed are the best-fit parameters for each of the four parametric models considered in this paper. The listed values are the median of 
the probability distribution function of each parameter. The errors correspond to the 16th and 84th percentiles around the median. The 
columns are: (1) parametric model. Models are listed in order of decreasing Bayesian evidence, (2) X-ray luminosity function normalisation 
(see equation]^, (3) break luminosity of the X-ray luminosity function (see equation]^, (4) faint-end slope, (5) bright-end slope, (6) 
density evolution parameter (see equations [^ [^ [^ , (7) luminosity evolution parameters (see equations , (8) /3 parameter for the 
LDDE model, (9) base 10 logarithm of the Bayesian evidence for each model, (10) the difference between the logj^Q Z of each model and 
the PDE that has the highest evidence. The bayes factor of the PDE model relative to each of the other three is exp(Aln.E). 


the redshift intervals 2 = 3 — 4 and 4 — 5. The analysis 
methodology we develop takes into account Poisson errors 
in the determination of X-ray fluxes and luminosities as well 
as uncertainties in photometric redshift measurements. We 
demonstrate that the latter is critical for unbiased measure¬ 
ments of the AGN space density at high redshift. Ignoring 
photometric redshift errors overestimates the X-ray lumi¬ 
nosity function. We also choose to follow a non-parametric 
approach and determine the space density of AGN in lumi¬ 
nosity and redshift bins. This allows us to explore the evolu¬ 
tion of the AGN population independent of model assump¬ 
tions. Additionally when a model parametrisation is applied 
to AGN at all redshifts the fit may be driven by the redshift 
and luminosity intervals that contain most data. This can 
introduce biases at high redshift, 2^3, where AGN samples 
are typically small as a result of the rapid evolution of the 
AGN population as well as survey sensitivity and volume 
limitations. We account for this potential issue by splitting 
the luminosity function into independent redshift compo¬ 
nents, i.e. 2 < 3, 2 = 3 — 5. 


We confirm previous studies for a strong evolution of 
the AGN population in the redshift interval 2 = 3 — 5 
rusa et al.||2009| |Vito et al.|[^13[ [Givano et al.|[2012[ 


Kalfountzou et al. 20141. We also find suggestive evidence 


for luminosity dependent evolution of the X-ray luminos¬ 
ity function. The space density of AGN in the luminos¬ 
ity interval log Lx (2 — lOkeV) « 43 — 45ergs“^ decreases 
faster than more luminous sources between redshifts 2 = 3.5 
and 2 = 4.5, albeit at the 90% significance level. A similar 
evolution pattern is also observed in the optical luminos¬ 
ity function of QSOs between redshifts 2 = 3 and 2 = 4 
( [Masters et al.|2012[ ). Our finding can be interpreted as evi¬ 
dence that the formation epoch of the most powerful QSOs 
[Lx(2 — lOkeV) > 10^® ergs“^] precedes that of lower lumi¬ 
nosity systems. This is similar to AGN downsizing trends 


established at lower redshifts (Ueda et al. 2003 Hasinger 
2008 Aird et ^[2010 Miyaji et al.||20f5 l. A strong evolu¬ 


tion of the faint-end of the AGN luminosity function with 
increasing redshift is consistent with the absence of X-ray se- 
lected AGN at a 5 in the GANDELS (Grogin et al.||2011 
Koekemoer et al.j^Oll I subregion of the (Jhandra Deep Eiefc 


South (Weigel et al.|20l5|. Extrapolating the LDDE model 


of Table we predict < 1 AGN at redshift 2 > 5 in that 
field. 


Our analysis also places limits on the contribution of 
AGN to the UV photon field needed to keep the hydro¬ 
gen ionised at high redshift. Using empirical relations for 
the type-1 AGN fraction as a function of luminosity ( jMer-j 
loni et al. 2014) we show that AGN dominate or at least 


contribute a sizable fraction of the required UV photons 
to redshift 2 « 4. At higher redshift the evolution of the 
X-ray luminosity function translates to a decreasing contri¬ 
bution of X-ray AGN to the UV photon field required to 
keep the hydrogen ionised. The extreme assumption of a 
photon escaping fraction of unity for all AGN sets an upper 
limit of 30% to the contribution of AGN to the UV photon 
rate density required to keep the hydrogen ionised at 2 = 5. 


Barger et al. (2003a I use multicolour optical data in the 2 Ms 


Chandra Deep Field North (Barger et al. 2003bI and con¬ 


clude that the X-ray selected AGN candidates at 2 = 5 — 6.5 
are too few to ionise the intergalactic medium at those red¬ 


shifts. Haardt & Madau (2012) estimated the contribution 


of AGN to hydrogen ionisation rate using the jUeda et al.j 
(20031 X-ray luminosity function and AGN obscuration dis¬ 


tribution. They found that AGN do not play an important 
role as a source of ionising photons above redshifts ~ 4. 


Grissom et al. (2014) determine the contribution of AGN to 


the ionisation of the hydrogen in the Universe by taking into 
account secondary collisional ionisations from the X-ray ra¬ 
diation. They extrapolate to high redshift (2 > 6) the 


Hiroi 


et al. (2012) hard X-ray luminosity function and conclude 


that AGN only contribute a small fraction of the photon 
rate densities required to ionise the Universe at these red¬ 
shifts. Our results are in agreement with the above studies 
and do not support claims for a dominant role of AGN to 
the ionisation of the hydrogen in the Universe at redshift 
2 > 4 (Glikman et al.|2011 Giallongo et al.|2015 1. This dis- 
crepancy is likely related to the way ditterent groups select 
their samples and subsequently account for this selection in 
the analysis. It also highlights the need for further research 
to better constrain the impact of AGN radiation to the ioni¬ 


sation of the Universe. Glikman et al. (20111 estimate type-I 


QSO space densities at 2 « 4 that are a factor of 3-4 higher 


than those determined by Masters et al. (20121 or Ikeda et al. 


(2011) at similar redshifts and luminosities. Giallongo et al. 
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120151 combined X-ray and HST optical/near-IR data in the 
CANDELS GOODS-S region to identify optical sources with 
photometric or spectroscopic redshifts z > 4 and then study 
their X-ray properties following methods described by [Fiore] 


et al. (20121. Their approach allows them to identify faint 
AGN with X-ray luminosities as a low as Lx ~ 10^^ ergs“^ 
in the redshift interval z = 4 — 6.5. Strictly speaking thejGi- 


allongo et al. (20151 photon-rate densities in Figure 14 are 


upper limits. The UV photon escape fraction is set to one, 
part of the observed UV radiation of individual sources may 
be associated with the AGN host galaxy, the photometric 
redshift uncertainties of at least some sources in the sample 
are large. The increase of the X-ray depth in the GANDELS 
GOODS-S from 4 to 7 Ms (PI Brandt) region will help better 
constrain the faint-end of the AGN luminosity function at 
high redshift and their role in the ionisation of the Universe. 

Finally the parametric X-ray luminosity functions de¬ 
rived in this paper are used to make predictions on the num¬ 
ber of z > 3 AGN that eROSITA ( [Merloni et al.|[2012[ ) sur¬ 
veys will detect. Our LDDE parametrisation is convolved 
with the expected X-ray sensitivity of the 4-year eROSITA 
All Sky Survey. The number of AGN in logarithmic lumi¬ 
nosity bins of size A log Lx = 0.5 is plotted as a function of 
2-10 keV luminosity in Figure [T^ Predictions are presented 
for 3 redshift intervals, 2 = 3 —4, 4 —5 and 5 — 6. This shows 
that surveys by eROSITA will provide tight constraints on 
the evolution of bright AGN and will allow us to explore 
with high statistical significance the evidence for luminosity 
dependence of the AGN population at high redshift. This 
however, would also require a dedicate follow-up program 
to identify high redshift AGN among the eROSITA popula¬ 
tion. High multiplex spectroscopic facilities that are able to 
simultaneously observe large number of targets over a wide 
field of view are essential for eROSITA X-ray source follow¬ 
ups. The SDSS/BOSS spectrogra phs (jSmee et ar]|2013 1 at 
the Apache Point SDSS telescope ( Gunn et al.|2006 1 and in 
the future the ESO/4MOST facility (4-metre Multi-Object 
Spectroscopic Telescope, de Jong et al.|20T4 1 are well suited 
for follow-up observations of the eROSITA sky. 


• the faint-end slope of UV/optical QSO luminosity func¬ 
tions is steeper than that of the X-ray selected AGN sam¬ 
ples. This implies an increasing fraction of type-I AGN with 
decreasing X-ray luminosity at 2 > 3. 

• X-ray AGN may dominate or at least contribute sub¬ 
stantially to the UV photon rate density required to keep 
the Universe ionised to 2 = 4. At higher redshift the contri¬ 
bution of AGN to UV hydrogen ionising field decreases. 
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6 CONCLUSIONS 

X-ray data from Chandra deep surveys and the shallow/wide 
XMM-XXL sample are combined to explore the evolution of 
the X-ray luminosity function at high redshift, 2 = 3 — 5. 
Our analysis accounts for Poisson errors in the calculation 
of fluxes and luminosities as well as photometric redshift 
uncertainties. We also show that the latter point is crucial 
for unbiased AGN space density measurements. The sample 
used in the paper consists of nearly 340 sources with either 
photometric (212) or spectroscopic (128) redshift in the red¬ 
shift range 2 = 3 — 5. The luminosity baseline of the sample 
is Lx (2 — lOkeV) « 10"*^ — 10^® erg s^ at 2 > 3. Our main 
findings are 

• the AGN population evolves strongly between the red¬ 
shift intervals 2 = 3 — 4 and 2 = 4 — 5 

• there is evidence, significant at the 90% level, that the 
amplitude of the AGN evolution depends on X-ray luminos¬ 
ity. Sources with luminosities Lx (2 — lOkeV) < 10^® ergs“^ 
appear to evolve faster than brighter ones. 
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